## DADOS DE FAMA-FRENCH:
# 1. Importação dos Dados
ff_raw <- read_csv("25_Portfolios_5x5.csv", skip = 15)
# 2. Limpeza e Tratamento Mensal
ff_monthly <- ff_raw %>%
rename(Date = `...1`) %>%
filter(!is.na(Date), !str_detect(Date, "Copyright")) %>%
mutate(Date = ymd(paste0(Date, "01"))) %>%
mutate(across(-Date, as.numeric))
# Carregar e preparar os dados de CPI
cpi_raw <- read_csv("CPIAUCNS.csv")
# Cálculo do componente inflacionário
cpi_monthly <- cpi_raw %>%
rename(CPI = CPIAUCNS) %>%
mutate(Date = ymd(observation_date)) %>%
arrange(Date) %>%
mutate(Inflation_pct = (CPI / lag(CPI) - 1) * 100) %>%
select(Date, Inflation_pct)
# Adicionando o CPI ao dataframe
ff_monthly_real <- ff_monthly %>%
left_join(cpi_monthly, by = "Date") %>%
filter(!is.na(Inflation_pct))
# 3. Conversão para Trimestral (compondo os retornos com o componente inflacionário)
ff_quarterly_real <- ff_monthly_real %>%
mutate(YearQuarter = as.yearqtr(Date)) %>%
group_by(YearQuarter) %>%
summarise(
across(
.cols = -c(Date, Inflation_pct),
.fns = ~ ((prod(1 + .x / 100, na.rm = TRUE) / prod(1 + Inflation_pct / 100, na.rm = TRUE)) - 1) * 100
),comp_inflation_quarterly = (prod(1 + Inflation_pct / 100, na.rm = TRUE) - 1) * 100
) %>%
ungroup()
head(ff_quarterly_real)
## # A tibble: 6 × 27
## YearQuarter `SMALL LoBM` `ME1 BM2` `ME1 BM3` `ME1 BM4` `SMALL HiBM` `ME2 BM1`
## <yearqtr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1947 Q1 40.8 26.9 51.6 71.1 114. 41.3
## 2 1947 Q2 66.6 46.8 107. 120. 205. 78.4
## 3 1947 Q3 54.7 71.6 119. 217. 374. 112.
## 4 1947 Q4 32.8 92.5 148. 273. 429. 125.
## 5 1948 Q1 88.4 104. 162. 295. 499. 138.
## 6 1948 Q2 73.9 149. 178. 302. 511. 162.
## # ℹ 20 more variables: `ME2 BM2` <dbl>, `ME2 BM3` <dbl>, `ME2 BM4` <dbl>,
## # `ME2 BM5` <dbl>, `ME3 BM1` <dbl>, `ME3 BM2` <dbl>, `ME3 BM3` <dbl>,
## # `ME3 BM4` <dbl>, `ME3 BM5` <dbl>, `ME4 BM1` <dbl>, `ME4 BM2` <dbl>,
## # `ME4 BM3` <dbl>, `ME4 BM4` <dbl>, `ME4 BM5` <dbl>, `BIG LoBM` <dbl>,
## # `ME5 BM2` <dbl>, `ME5 BM3` <dbl>, `ME5 BM4` <dbl>, `BIG HiBM` <dbl>,
## # comp_inflation_quarterly <dbl>
## DADOS DE CONSUMO:
# 1. Importação do Consumo Nominal (Tabela 2.3.5)
# Carregar os cabeçalhos
header_nominal <- read_csv("Table 2.3.5. Personal Consumption Expenditures by Major Type of Product.csv",
skip = 3,
n_max = 2,
col_names = FALSE
)
# Carregar os dados
data_nominal <- read_csv(
"Table 2.3.5. Personal Consumption Expenditures by Major Type of Product.csv",
skip = 5,
col_names = FALSE,
col_types = cols(.default = "c")
)
# Criar nomes de colunas corretos
year_row_nom <- header_nominal[1, -c(1, 2)]
quarter_row_nom <- header_nominal[2, -c(1, 2)]
new_column_names_nom <- paste0(zoo::na.locf(unlist(year_row_nom)), "_", unlist(quarter_row_nom))
colnames(data_nominal) <- c("Line", "Description", new_column_names_nom)
# Isolar a série de consumo nominal de não duráveis (Linha 8)
nondurables_nominal <- data_nominal %>%
filter(Line == "8") %>%
pivot_longer(
cols = contains("_Q"),
names_to = "Quarter_str",
values_to = "Nondurables_Nominal"
) %>%
mutate(
YearQuarter = as.yearqtr(Quarter_str, format = "%Y_Q%q"),
Nondurables_Nominal = as.numeric(Nondurables_Nominal)
) %>%
select(YearQuarter, Nondurables_Nominal)
# 2. Importação do Índice de Preços (Tabela 2.3.4)
# Carregar os cabeçalhos
header_prices <- read_csv("Table 2.3.4. Price Indexes for Personal Consumption Expenditures by Major Type of Product.csv",
skip = 3,
n_max = 2,
col_names = FALSE
)
# Carregar os dados
data_prices <- read_csv(
"Table 2.3.4. Price Indexes for Personal Consumption Expenditures by Major Type of Product.csv",
skip = 5,
col_names = FALSE,
col_types = cols(.default = "c")
)
# Criar nomes de colunas corretos
year_row_price <- header_prices[1, -c(1, 2)]
quarter_row_price <- header_prices[2, -c(1, 2)]
new_column_names_price <- paste0(zoo::na.locf(unlist(year_row_price)), "_", unlist(quarter_row_price))
colnames(data_prices) <- c("Line", "Description", new_column_names_price)
# Isolar a série do índice de preços de não duráveis (Linha 8)
price_index_nondurables <- data_prices %>%
filter(Line == "8") %>%
pivot_longer(
cols = contains("_Q"),
names_to = "Quarter_str",
values_to = "Price_Index"
) %>%
mutate(
YearQuarter = as.yearqtr(Quarter_str, format = "%Y_Q%q"),
Price_Index = as.numeric(Price_Index)
) %>%
select(YearQuarter, Price_Index)
# 3. Calculando o Consumo Real
# Unir as duas tabelas e deflacionar os valores nominais
nondurables_quarterly <- nondurables_nominal %>%
inner_join(price_index_nondurables, by = "YearQuarter") %>%
mutate(
# A fórmula para deflacionar é (Valor Nominal / Índice de Preços) * 100
Nondurables_Real = (Nondurables_Nominal / Price_Index) * 100
) %>%
select(YearQuarter, Nondurables_Real)
cat("--- Amostra do Consumo Real de Não-Duráveis Calculado ---\n")
## --- Amostra do Consumo Real de Não-Duráveis Calculado ---
head(nondurables_quarterly)
## # A tibble: 6 × 2
## YearQuarter Nondurables_Real
## <yearqtr> <dbl>
## 1 1947 Q1 502.
## 2 1947 Q2 510.
## 3 1947 Q3 512.
## 4 1947 Q4 506.
## 5 1948 Q1 507.
## 6 1948 Q2 513.
# 4. Importação e Adequação do Dataframe de População (Tabela 2.1)
header_pop <- read_csv("Table 2.1. Personal Income and Its Disposition.csv",
skip = 3,
n_max = 2,
col_names = FALSE
)
data_pop <- read_csv(
"Table 2.1. Personal Income and Its Disposition.csv",
skip = 5,
col_names = FALSE,
col_types = cols(.default = "c")
)
# Criando os nomes de coluna corretos
year_row_pop <- header_pop[1, -c(1, 2)]
quarter_row_pop <- header_pop[2, -c(1, 2)]
new_column_names_pop <- paste0(zoo::na.locf(unlist(year_row_pop)), "_", unlist(quarter_row_pop))
colnames(data_pop) <- c("Line", "Description", new_column_names_pop)
# 5. Limpeza e Transformação para a série de População
population_quarterly <- data_pop %>%
filter(Line == 40) %>%
pivot_longer(
cols = contains("_Q"),
names_to = "Quarter_str",
values_to = "Population_Thousands"
) %>%
mutate(
YearQuarter = as.yearqtr(Quarter_str, format = "%Y_Q%q"),
Population_Thousands = as.numeric(Population_Thousands)
) %>%
select(YearQuarter, Population_Thousands)
# 6. Unindo os dataframes de consumo e população
consumption_full <- nondurables_quarterly %>%
left_join(population_quarterly, by = "YearQuarter")
# Verificando se há NAs na nova coluna
cat("\nValores Faltantes (NA) na coluna de População:", sum(is.na(consumption_full$Population_Thousands)), "\n")
##
## Valores Faltantes (NA) na coluna de População: 0
# 7. Criando a coluna de consumo per-capita
consumption_full <- consumption_full %>%
mutate(Nondurables_Real_PC = (Nondurables_Real * 1000000000) / (Population_Thousands * 1000))
summary(consumption_full$Nondurables_Real_PC)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3446 4496 5850 6030 7686 10256
head(consumption_full)
## # A tibble: 6 × 4
## YearQuarter Nondurables_Real Population_Thousands Nondurables_Real_PC
## <yearqtr> <dbl> <dbl> <dbl>
## 1 1947 Q1 502. 143143 3505.
## 2 1947 Q2 510. 143790 3548.
## 3 1947 Q3 512. 144449 3546.
## 4 1947 Q4 506. 145122 3487.
## 5 1948 Q1 507. 145709 3481.
## 6 1948 Q2 513. 146289 3504.
# 8. Calculando o Crescimento do Consumo Real Per-Capita
consumption_final <- consumption_full %>%
arrange(YearQuarter) %>%
mutate(
consumption_growth_pc = log(Nondurables_Real_PC) - lag(log(Nondurables_Real_PC))
)
head(consumption_final)
## # A tibble: 6 × 5
## YearQuarter Nondurables_Real Population_Thousands Nondurables_Real_PC
## <yearqtr> <dbl> <dbl> <dbl>
## 1 1947 Q1 502. 143143 3505.
## 2 1947 Q2 510. 143790 3548.
## 3 1947 Q3 512. 144449 3546.
## 4 1947 Q4 506. 145122 3487.
## 5 1948 Q1 507. 145709 3481.
## 6 1948 Q2 513. 146289 3504.
## # ℹ 1 more variable: consumption_growth_pc <dbl>
## DADOS DOS FATORES FAMA-FRENCH E TAXA LIVRE DE RISCO (RF):
# 1. Importação dos Dados
factors_raw <- read_csv("F-F_Research_Data_Factors.csv", skip = 3)
# 2. Limpeza e Tratamento Mensal
factors_monthly <- factors_raw %>%
rename(Date = `...1`) %>%
filter(!is.na(Date)) %>%
mutate(Date = ymd(paste0(Date, "01"))) %>%
mutate(across(-Date, as.numeric))
# 3. Conversão da RF para Trimestral (compondo a taxa)
rf_quarterly <- factors_monthly %>%
mutate(YearQuarter = as.yearqtr(Date)) %>%
group_by(YearQuarter) %>%
summarise(
RF_quarterly = (prod(1 + RF / 100, na.rm = TRUE) - 1) * 100
) %>%
ungroup()
## JUNÇÃO FINAL E CÁLCULO DO SDF
# 1. Unir as bases de dados (Fama-French, Consumo e Taxa Livre de Risco)
final_data <- ff_quarterly_real %>%
left_join(consumption_final, by = "YearQuarter") %>%
left_join(rf_quarterly, by = "YearQuarter")
# 2. Criar o Dataframe Final de Análise
analysis_df <- final_data %>%
# 2.1 Calcular o excesso de retorno (prêmio de risco) para cada portfólio
mutate(across(
.cols = starts_with("SMALL") | starts_with("ME") | starts_with("BIG"),
.fns = ~ .x - RF_quarterly,
.names = "{.col}_excess"
)) %>%
# 2.2 Calcular o SDF
mutate(
gross_consumption_growth = exp(consumption_growth_pc),
SDF = 0.99 * (gross_consumption_growth ^ -2)
) %>%
# 2.3 Remover a primeira linha que tem NA no SDF (devido ao lag do crescimento)
filter(!is.na(SDF))
# 3. Calcular o Prêmio de Risco Estimado pelo Modelo
mean_sdf <- mean(analysis_df$SDF, na.rm = TRUE)
model_premiums <- analysis_df %>%
summarise(across(
.cols = ends_with("_excess"),
.fns = ~ -cov(SDF, .x, use = "complete.obs") / mean_sdf,
.names = "model_{.col}"
))
model_premiums_long <- model_premiums %>%
pivot_longer(cols = everything(),
names_to = "Portfolio",
values_to = "Model_Estimated_Premium")
cat("--- Prêmios de Risco Trimestrais Estimados pelo Modelo\n")
## --- Prêmios de Risco Trimestrais Estimados pelo Modelo
print(model_premiums_long, n = 25)
## # A tibble: 25 × 2
## Portfolio Model_Estimated_Premium
## <chr> <dbl>
## 1 model_SMALL LoBM_excess 490.
## 2 model_SMALL HiBM_excess 497.
## 3 model_ME1 BM2_excess 2.63
## 4 model_ME1 BM3_excess 63.1
## 5 model_ME1 BM4_excess 182.
## 6 model_ME2 BM1_excess 1177.
## 7 model_ME2 BM2_excess 2260.
## 8 model_ME2 BM3_excess 631.
## 9 model_ME2 BM4_excess 2110.
## 10 model_ME2 BM5_excess 421.
## 11 model_ME3 BM1_excess 11420.
## 12 model_ME3 BM2_excess 23340.
## 13 model_ME3 BM3_excess 14701.
## 14 model_ME3 BM4_excess 16046.
## 15 model_ME3 BM5_excess 3249.
## 16 model_ME4 BM1_excess 414332.
## 17 model_ME4 BM2_excess 222608.
## 18 model_ME4 BM3_excess 145208.
## 19 model_ME4 BM4_excess 122116.
## 20 model_ME4 BM5_excess 66316.
## 21 model_ME5 BM2_excess 154703655.
## 22 model_ME5 BM3_excess 17136945.
## 23 model_ME5 BM4_excess 30761523.
## 24 model_BIG LoBM_excess 1210029838.
## 25 model_BIG HiBM_excess -8191002.
# 1. Tratando os Dados dos Fatores (análogo ao tratamento para os 25 portifólios da letra (a))
factors_monthly_raw <- read_csv("F-F_Research_Data_Factors.csv", skip = 3)
factors_monthly <- factors_monthly_raw %>%
rename(Date = `...1`) %>%
filter(!is.na(Date), !str_detect(Date, "Copyright")) %>%
mutate(Date = ymd(paste0(Date, "01"))) %>%
mutate(across(-Date, as.numeric))
factors_monthly_real <- factors_monthly %>%
left_join(cpi_monthly, by = "Date") %>%
filter(!is.na(Inflation_pct))
factors_quarterly_real <- factors_monthly_real %>%
mutate(YearQuarter = as.yearqtr(Date)) %>%
group_by(YearQuarter) %>%
summarise(
across(
.cols = c(`Mkt-RF`, RF),
.fns = ~ ((prod(1 + .x / 100, na.rm = TRUE) / prod(1 + Inflation_pct / 100, na.rm = TRUE)) - 1) * 100,
.names = "{.col}_real"
),
`Mkt-RF` = (prod(1 + `Mkt-RF` / 100, na.rm = TRUE) - 1) * 100,
SMB = (prod(1 + SMB / 100, na.rm = TRUE) - 1) * 100,
HML = (prod(1 + HML / 100, na.rm = TRUE) - 1) * 100
) %>%
ungroup()
cat("--- Fatores Fama-French Trimestrais (Reais e Nominais) ---\n")
## --- Fatores Fama-French Trimestrais (Reais e Nominais) ---
head(factors_quarterly_real)
## # A tibble: 6 × 6
## YearQuarter `Mkt-RF_real` RF_real `Mkt-RF` SMB HML
## <yearqtr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1947 Q1 -4.53 -1.77 -2.75 -0.911 0.751
## 2 1947 Q2 -1.17 -0.365 -0.717 -7.42 0.525
## 3 1947 Q3 -2.66 -4.23 1.77 3.45 4.30
## 4 1947 Q4 1.69 -1.51 3.45 -3.62 4.97
## 5 1948 Q1 -0.688 0.230 -0.688 0.898 6.02
## 6 1948 Q2 7.89 -2.66 11.1 -2.55 5.72
# 2. Unindo os Retornos dos Portfólios com o Fator de Mercado
capm_data <- analysis_df %>%
select(YearQuarter, ends_with("_excess")) %>%
inner_join(factors_quarterly_real %>% select(YearQuarter, `Mkt-RF`), by = "YearQuarter")
# 3. Estimando o CAPM para cada um dos 25 Portfólios
capm_results <- capm_data %>%
pivot_longer(
cols = ends_with("_excess"),
names_to = "Portfolio",
values_to = "Excess_Return"
) %>%
# Agrupa por portfólio para rodar uma regressão para cada um
group_by(Portfolio) %>%
# tidyr::nest() cria um "sub-dataframe" para cada portfólio
nest() %>%
# purrr::map() aplica a função de regressão (lm) a cada sub-dataframe
mutate(
model = map(data, ~ lm(Excess_Return ~ `Mkt-RF`, data = .x)),
# broom::tidy() extrai os coeficientes (alpha e beta) de cada modelo de forma organizada
tidied = map(model, tidy)
) %>%
# Expande a lista de coeficientes de volta para o formato de tabela
unnest(tidied)
cat("\n\n--- Coeficientes do CAPM (Alpha e Beta) para cada Portfólio ---\n")
##
##
## --- Coeficientes do CAPM (Alpha e Beta) para cada Portfólio ---
print(capm_results)
## # A tibble: 50 × 8
## # Groups: Portfolio [25]
## Portfolio data model term estimate std.error statistic p.value
## <chr> <list> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 SMALL LoBM_excess <tibble> <lm> (Int… 127810. 12113. 10.6 1.88e-22
## 2 SMALL LoBM_excess <tibble> <lm> `Mkt… 217. 1273. 0.170 8.65e- 1
## 3 SMALL HiBM_excess <tibble> <lm> (Int… 375921. 28416. 13.2 5.31e-32
## 4 SMALL HiBM_excess <tibble> <lm> `Mkt… 8194. 2986. 2.74 6.41e- 3
## 5 ME1 BM2_excess <tibble> <lm> (Int… 69114. 6924. 9.98 1.54e-20
## 6 ME1 BM2_excess <tibble> <lm> `Mkt… 1460. 727. 2.01 4.56e- 2
## 7 ME1 BM3_excess <tibble> <lm> (Int… 79808. 7416. 10.8 3.61e-23
## 8 ME1 BM3_excess <tibble> <lm> `Mkt… 1647. 779. 2.11 3.54e- 2
## 9 ME1 BM4_excess <tibble> <lm> (Int… 141324. 11955. 11.8 7.02e-27
## 10 ME1 BM4_excess <tibble> <lm> `Mkt… 2961. 1256. 2.36 1.90e- 2
## # ℹ 40 more rows
# 4. Calculando o Prêmio de Risco Previsto pelo CAPM
# Calculamos a média do prêmio de risco de mercado histórico
avg_market_premium <- mean(capm_data$`Mkt-RF`, na.rm = TRUE)
# Filtramos para pegar apenas os Betas (o coeficiente de `Mkt-RF`) de cada portfólio
betas <- capm_results %>%
filter(term == "`Mkt-RF`") %>%
select(Portfolio, beta = estimate)
# Calculamos o prêmio de risco de cada portfólio: Beta * Média do Prêmio de Mercado
capm_risk_premiums <- betas %>%
mutate(CAPM_Premium = beta * avg_market_premium)
cat("\n\n--- Prêmio de Risco Trimestral Estimado pelo CAPM\n")
##
##
## --- Prêmio de Risco Trimestral Estimado pelo CAPM
print(capm_risk_premiums, n = 25)
## # A tibble: 25 × 3
## # Groups: Portfolio [25]
## Portfolio beta CAPM_Premium
## <chr> <dbl> <dbl>
## 1 SMALL LoBM_excess 217. 487.
## 2 SMALL HiBM_excess 8194. 18390.
## 3 ME1 BM2_excess 1460. 3278.
## 4 ME1 BM3_excess 1647. 3696.
## 5 ME1 BM4_excess 2961. 6645.
## 6 ME2 BM1_excess 35473. 79616.
## 7 ME2 BM2_excess 46909. 105283.
## 8 ME2 BM3_excess 27898. 62613.
## 9 ME2 BM4_excess 30516. 68490.
## 10 ME2 BM5_excess 18057. 40526.
## 11 ME3 BM1_excess 477641. 1072009.
## 12 ME3 BM2_excess 598170. 1342521.
## 13 ME3 BM3_excess 236064. 529817.
## 14 ME3 BM4_excess 221389. 496881.
## 15 ME3 BM5_excess 77299. 173488.
## 16 ME4 BM1_excess 14184253. 31834846.
## 17 ME4 BM2_excess 4805607. 10785607.
## 18 ME4 BM3_excess 2732268. 6132247.
## 19 ME4 BM4_excess 2320861. 5208892.
## 20 ME4 BM5_excess 1379408. 3095915.
## 21 ME5 BM2_excess 2616557520. 5872548148.
## 22 ME5 BM3_excess 359205227. 806192860.
## 23 ME5 BM4_excess 986420947. 2213902985.
## 24 BIG LoBM_excess 17557041403. 39404664404.
## 25 BIG HiBM_excess 214126499. 480581132.
# 1. Preparando os Dados para a Regressão de 3 Fatores
# Unir os retornos reais dos portfólios com os três fatores Fama-French.
three_factor_data <- analysis_df %>%
select(YearQuarter, ends_with("_excess")) %>%
inner_join(
factors_quarterly_real %>% select(YearQuarter, `Mkt-RF_real`, SMB, HML),
by = "YearQuarter"
)
cat("--- Amostra dos Dados para o Modelo de 3 Fatores ---\n")
## --- Amostra dos Dados para o Modelo de 3 Fatores ---
head(three_factor_data)
## # A tibble: 6 × 29
## YearQuarter `SMALL LoBM_excess` `SMALL HiBM_excess` `ME1 BM2_excess`
## <yearqtr> <dbl> <dbl> <dbl>
## 1 1947 Q2 66.5 205. 46.7
## 2 1947 Q3 54.6 374. 71.5
## 3 1947 Q4 32.6 428. 92.3
## 4 1948 Q1 88.2 499. 104.
## 5 1948 Q2 73.7 511. 148.
## 6 1948 Q3 2.34 320. 50.5
## # ℹ 25 more variables: `ME1 BM3_excess` <dbl>, `ME1 BM4_excess` <dbl>,
## # `ME2 BM1_excess` <dbl>, `ME2 BM2_excess` <dbl>, `ME2 BM3_excess` <dbl>,
## # `ME2 BM4_excess` <dbl>, `ME2 BM5_excess` <dbl>, `ME3 BM1_excess` <dbl>,
## # `ME3 BM2_excess` <dbl>, `ME3 BM3_excess` <dbl>, `ME3 BM4_excess` <dbl>,
## # `ME3 BM5_excess` <dbl>, `ME4 BM1_excess` <dbl>, `ME4 BM2_excess` <dbl>,
## # `ME4 BM3_excess` <dbl>, `ME4 BM4_excess` <dbl>, `ME4 BM5_excess` <dbl>,
## # `ME5 BM2_excess` <dbl>, `ME5 BM3_excess` <dbl>, `ME5 BM4_excess` <dbl>, …
# 2. Estimando o Modelo de 3 Fatores para cada Portfólio (análogo à letra (b))
three_factor_results <- three_factor_data %>%
pivot_longer(
cols = ends_with("_excess"),
names_to = "Portfolio",
values_to = "Excess_Return"
) %>%
group_by(Portfolio) %>%
nest() %>%
mutate(
model = map(data, ~ lm(Excess_Return ~ `Mkt-RF_real` + SMB + HML, data = .x)),
tidied = map(model, tidy)
) %>%
unnest(tidied)
cat("\n\n--- Coeficientes do Modelo de 3 Fatores para cada Portfólio ---\n")
##
##
## --- Coeficientes do Modelo de 3 Fatores para cada Portfólio ---
print(three_factor_results, n = 12)
## # A tibble: 100 × 8
## # Groups: Portfolio [25]
## Portfolio data model term estimate std.error statistic p.value
## <chr> <list> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 SMALL LoBM_excess <tibble> <lm> (Int… 133129. 12060. 11.0 4.17e-24
## 2 SMALL LoBM_excess <tibble> <lm> `Mkt… 154. 1315. 0.117 9.07e- 1
## 3 SMALL LoBM_excess <tibble> <lm> SMB -1082. 2173. -0.498 6.19e- 1
## 4 SMALL LoBM_excess <tibble> <lm> HML -4456. 2013. -2.21 2.76e- 2
## 5 SMALL HiBM_excess <tibble> <lm> (Int… 378767. 28449. 13.3 2.82e-32
## 6 SMALL HiBM_excess <tibble> <lm> `Mkt… 8756. 3103. 2.82 5.08e- 3
## 7 SMALL HiBM_excess <tibble> <lm> SMB 1112. 5127. 0.217 8.28e- 1
## 8 SMALL HiBM_excess <tibble> <lm> HML 2867. 4749. 0.604 5.47e- 1
## 9 ME1 BM2_excess <tibble> <lm> (Int… 70548. 6943. 10.2 4.09e-21
## 10 ME1 BM2_excess <tibble> <lm> `Mkt… 1404. 757. 1.85 6.47e- 2
## 11 ME1 BM2_excess <tibble> <lm> SMB 500. 1251. 0.400 6.90e- 1
## 12 ME1 BM2_excess <tibble> <lm> HML -318. 1159. -0.274 7.84e- 1
## # ℹ 88 more rows
# 3. Calculando o Prêmio de Risco Previsto pelo Modelo de 3 Fatores
# Calculamos a média histórica de cada um dos três fatores
avg_factor_premiums <- three_factor_data %>%
summarise(
avg_mkt_real = mean(`Mkt-RF_real`, na.rm = TRUE),
avg_smb = mean(SMB, na.rm = TRUE),
avg_hml = mean(HML, na.rm = TRUE)
)
# Extraímos os betas de cada portfólio para cada fator
betas_3f <- three_factor_results %>%
filter(term %in% c("`Mkt-RF_real`", "SMB", "HML")) %>%
select(Portfolio, term, estimate) %>%
pivot_wider(names_from = term, values_from = estimate) %>%
rename(
beta_mkt = `\`Mkt-RF_real\``,
beta_smb = SMB,
beta_hml = HML
)
# Calculamos o prêmio de risco para cada portfólio
three_factor_premiums <- betas_3f %>%
mutate(
Three_Factor_Premium =
beta_mkt * avg_factor_premiums$avg_mkt_real +
beta_smb * avg_factor_premiums$avg_smb +
beta_hml * avg_factor_premiums$avg_hml
)
cat("\n\n--- Prêmio de Risco Trimestral Estimado pelo Modelo de 3 Fatores\n")
##
##
## --- Prêmio de Risco Trimestral Estimado pelo Modelo de 3 Fatores
print(three_factor_premiums, n = 25)
## # A tibble: 25 × 5
## # Groups: Portfolio [25]
## Portfolio beta_mkt beta_smb beta_hml Three_Factor_Premium
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 SMALL LoBM_excess 154. -1082. -4.46e3 -4832.
## 2 SMALL HiBM_excess 8756. 1112. 2.87e3 15544.
## 3 ME1 BM2_excess 1404. 500. -3.18e2 1844.
## 4 ME1 BM3_excess 1619. 386. -7.22e2 1677.
## 5 ME1 BM4_excess 2073. 6802. 4.13e2 6388.
## 6 ME2 BM1_excess 36715. -4253. -7.11e3 41571.
## 7 ME2 BM2_excess 45972. 4310. 2.59e3 68188.
## 8 ME2 BM3_excess 30694. -12397. 2.63e3 39490.
## 9 ME2 BM4_excess 31687. -3243. 7.59e3 50092.
## 10 ME2 BM5_excess 20864. -11356. 1.29e4 36873.
## 11 ME3 BM1_excess 484325. -31923. 6.30e4 719604.
## 12 ME3 BM2_excess 587048. 65336. 1.50e5 995153.
## 13 ME3 BM3_excess 243748. -51657. 3.23e4 346590.
## 14 ME3 BM4_excess 224658. -2379. 3.49e4 345264.
## 15 ME3 BM5_excess 94616. -85674. 5.66e4 149613.
## 16 ME4 BM1_excess 13808157. 1505782. 1.59e6 21410852.
## 17 ME4 BM2_excess 5002452. -1267525. 6.04e5 6957196.
## 18 ME4 BM3_excess 2850086. -340073. 1.60e6 5416120.
## 19 ME4 BM4_excess 2565810. -1424794. 6.56e5 3568838.
## 20 ME4 BM5_excess 1632878. -1239227. 1.56e6 3285212.
## 21 ME5 BM2_excess 2665886771. -387086610. 8.43e8 4371596994.
## 22 ME5 BM3_excess 456205123. -541643157. -1.75e8 205289163.
## 23 ME5 BM4_excess 1104354255. -633080922. 3.57e8 1603921830.
## 24 BIG LoBM_excess 20678555373. -20940905374. -4.96e9 13990816776.
## 25 BIG HiBM_excess 316906469. -546665652. 1.39e8 330948133.
ipca<-ipeadata("PRECOS12_IPCA12")%>%filter(date>=as.Date("2000-01-01"),date<=as.Date("2025-01-01"))%>%arrange(date)
saveRDS(ipca,"ipca_precos12_2000_2025.rds");write.csv(ipca,"ipca_precos12_2000_2025.csv",row.names=FALSE)
ipca<-readRDS("ipca_precos12_2000_2025.rds");s<-ts(ipca$value,start=c(2000,1),frequency=12)
summary(ipca$value);c(mean=mean(ipca$value),sd=sd(ipca$value),min=min(ipca$value),max=max(ipca$value))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1598 2574 3498 3858 5093 7112
## mean sd min max
## 3858.085 1577.834 1598.410 7111.860
ggplotly(autoplot(s)+theme_classic()+ggtitle("IPCA - Índice (2000-2025)"))%>%layout(plot_bgcolor="white",paper_bgcolor="white",yaxis=list(rangemode="tozero"))
ggplotly(autoplot(log(s))+theme_classic()+ggtitle("Log do IPCA"))%>%layout(plot_bgcolor="white",paper_bgcolor="white",yaxis=list(rangemode="tozero"))
ggplotly(autoplot(diff(log(s)))+theme_classic()+ggtitle("Δlog(IPCA)"))%>%layout(plot_bgcolor="white",paper_bgcolor="white",yaxis=list(rangemode="tozero"))
Pode-se observar uma notável tendência de crescimento do IPCA tanto em nível quanto em log, enquanto o comportamento dos dados em diferença do log apresentam um comportamento mais oscilatório, sem tendência clara de crescimento e se aproximando um pouco mais de um ruído branco.
adf.test(s);pp.test(s)
##
## Augmented Dickey-Fuller Test
##
## data: s
## Dickey-Fuller = -0.43059, Lag order = 6, p-value = 0.9846
## alternative hypothesis: stationary
##
## Phillips-Perron Unit Root Test
##
## data: s
## Dickey-Fuller Z(alpha) = -0.29147, Truncation lag parameter = 5,
## p-value = 0.99
## alternative hypothesis: stationary
Pode-se observar que os testes não rejeitam a hipótese nula de raiz unitária, dado o elevado p-valor dos testes, o que indica uma iminente tendência estocástica e não-estacionariedade, conforme indicado pelos gráficos.
adf.test(log(s));pp.test(log(s))
##
## Augmented Dickey-Fuller Test
##
## data: log(s)
## Dickey-Fuller = -2.3204, Lag order = 6, p-value = 0.4413
## alternative hypothesis: stationary
##
## Phillips-Perron Unit Root Test
##
## data: log(s)
## Dickey-Fuller Z(alpha) = -7.4022, Truncation lag parameter = 5, p-value
## = 0.6953
## alternative hypothesis: stationary
Assim como em nível, em log se observa um p-valor elevado, com os testes não rejeitando a hipótese nula de raiz unitária, indicando novamente tendência estocástica e não-estacionariedade
adf.test(diff(log(s)));pp.test(diff(log(s)))
##
## Augmented Dickey-Fuller Test
##
## data: diff(log(s))
## Dickey-Fuller = -5.813, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
##
## Phillips-Perron Unit Root Test
##
## data: diff(log(s))
## Dickey-Fuller Z(alpha) = -118, Truncation lag parameter = 5, p-value =
## 0.01
## alternative hypothesis: stationary
Para os dados em diferença do log, observa-se um p-valor consideravelmente menor, rejeitando a hipótese nula de raiz unitária e indicando estacionariedade.
d<-diff(log(s));a<-acf(d,plot=FALSE,na.action=na.pass);p<-pacf(d,plot=FALSE,na.action=na.pass);ci<-qnorm(0.975)/sqrt(a$n.used)
xa<-as.numeric(a$lag)*12;ya<-as.numeric(a$acf);xp<-as.numeric(p$lag)*12;yp<-as.numeric(p$acf)
plot_ly(x=xa,y=ya,type="bar")%>%layout(title="ACF Δlog(IPCA)",plot_bgcolor="white",paper_bgcolor="white",xaxis=list(title="Defasagem (meses)"),yaxis=list(title="Correlação",rangemode="tozero"),shapes=list(list(type="line",xref="x",yref="y",x0=min(xa),x1=max(xa),y0=ci,y1=ci),list(type="line",xref="x",yref="y",x0=min(xa),x1=max(xa),y0=-ci,y1=-ci)))
plot_ly(x=xp,y=yp,type="bar")%>%layout(title="PACF Δlog(IPCA)",plot_bgcolor="white",paper_bgcolor="white",xaxis=list(title="Defasagem (meses)"),yaxis=list(title="Correlação",rangemode="tozero"),shapes=list(list(type="line",xref="x",yref="y",x0=min(xp),x1=max(xp),y0=ci,y1=ci),list(type="line",xref="x",yref="y",x0=min(xp),x1=max(xp),y0=-ci,y1=-ci)))
Enquanto na PACF não se observa nenhum lag significativo, com todos os lags após o lag nulo se mantendo dentro do intervalo de confiança, pode-se assumir que as partes de média móvel sejam 0 tanto no aspecto não-sazonal quanto no sazonal, com possivelmente \(q=Q=0\).
Contudo, analisando a ACF, é possível observar um decaimento gradual após o lag 0 e além disso um comportamento de crescimento após um ano, indicando um possível fator AR tanto sazonal quanto não-sazonal, indicando possivelmente que \(p=P=1\).
Dessa forma, o modelo que tende a ser o mais adequado para o IPCA é o \(SARIMA(1,1,0)(1,0,0)_{12}\).
library(forecast);library(plotly)
y<-log(s)
fit1<-Arima(y,order=c(1,1,0),seasonal=list(order=c(1,0,0),period=12),include.drift=FALSE,method="ML")
fit2<-Arima(y,order=c(1,1,0),seasonal=list(order=c(0,0,1),period=12),include.drift=FALSE,method="ML")
summary(fit1);summary(fit2)
## Series: y
## ARIMA(1,1,0)(1,0,0)[12]
##
## Coefficients:
## ar1 sar1
## 0.8060 0.2269
## s.e. 0.0375 0.0656
##
## sigma^2 = 1.06e-05: log likelihood = 1295.33
## AIC=-2584.65 AICc=-2584.57 BIC=-2573.54
##
## Training set error measures:
## ME RMSE MAE MPE MAPE
## Training set 0.0007735025 0.003238879 0.002341851 0.009619096 0.02868748
## MASE ACF1
## Training set 0.03888154 -0.1301465
## Series: y
## ARIMA(1,1,0)(0,0,1)[12]
##
## Coefficients:
## ar1 sma1
## 0.8232 0.1775
## s.e. 0.0336 0.0559
##
## sigma^2 = 1.069e-05: log likelihood = 1294.07
## AIC=-2582.14 AICc=-2582.05 BIC=-2571.02
##
## Training set error measures:
## ME RMSE MAE MPE MAPE
## Training set 0.0007725472 0.003253154 0.002351729 0.009602986 0.02882324
## MASE ACF1
## Training set 0.03904554 -0.138668
Estimados os dois modelos, podemos observar valores muito próximos nos testes de AIC e BIC, mas o modelo com termo AR sazonal apresenta um valor substancialmente menor em ambos os testes, sendo esse o modelo mais indicado com base nos dois testes.
plres<-function(f,t){
r<-na.omit(residuals(f));ix<-seq_along(r)
a<-acf(r,plot=FALSE);p<-pacf(r,plot=FALSE);ci<-qnorm(0.975)/sqrt(a$n.used)
xa<-as.numeric(a$lag)*12;ya<-as.numeric(a$acf);xp<-as.numeric(p$lag)*12;yp<-as.numeric(p$acf)
rt<-plot_ly(x=ix,y=as.numeric(r),type="scatter",mode="lines")%>%layout(title=paste0("Resíduos - ",t),plot_bgcolor="white",paper_bgcolor="white",xaxis=list(title="Índice temporal"),yaxis=list(title="Resíduo",rangemode="tozero"))
ac<-plot_ly(x=xa,y=ya,type="bar")%>%layout(title=paste0("ACF dos resíduos - ",t),plot_bgcolor="white",paper_bgcolor="white",xaxis=list(title="Defasagem (meses)"),yaxis=list(title="Correlação",rangemode="tozero"),shapes=list(list(type="line",xref="x",yref="y",x0=min(xa),x1=max(xa),y0=ci,y1=ci),list(type="line",xref="x",yref="y",x0=min(xa),x1=max(xa),y0=-ci,y1=-ci)))
pc<-plot_ly(x=xp,y=yp,type="bar")%>%layout(title=paste0("PACF dos resíduos - ",t),plot_bgcolor="white",paper_bgcolor="white",xaxis=list(title="Defasagem (meses)"),yaxis=list(title="Correlação",rangemode="tozero"),shapes=list(list(type="line",xref="x",yref="y",x0=min(xp),x1=max(xp),y0=ci,y1=ci),list(type="line",xref="x",yref="y",x0=min(xp),x1=max(xp),y0=-ci,y1=-ci)))
list(rt=rt,ac=ac,pc=pc)
}
g1<-plres(fit1,"SARIMA(1,1,0)(1,0,0)[12]")
g2<-plres(fit2,"SARIMA(1,1,0)(0,0,1)[12]")
g1$rt;g2$rt
g1$ac;g2$ac
g1$pc;g2$pc
Dado que os dois modelos tem coeficientes próximos, é esperado um desempenho similar dos gráficos de resíduos, algo que de fato ocorre, com ambos aproximando-se de um ruído branco. Em ambos os modelos se observa um decaimento após o lag 0 com os lags \(1\) e \(2\) sendo negativos e fora do intervalo de confiança por uma diferença mínima além de um lag \(15\) significativo, o que poderia talvez indicar a necessidade de estimação de mais parâmetros para um ajuste mais fino, mas que perderia em parcimônia com a estimação de novos parâmetros. O mesmo vale para os valores dos lags \(1\) e \(2\) fora do intervalo de confiança da PACF em ambos os modelos.
lb<-function(f,k)c(LB24=Box.test(residuals(f),type="Ljung-Box",lag=24,fitdf=k)$p.value,
LB36=Box.test(residuals(f),type="Ljung-Box",lag=36,fitdf=k)$p.value)
rbind("(1,1,0)(1,0,0)[12]"=lb(fit1,2),
"(1,1,0)(0,0,1)[12]"=lb(fit2,2))
## LB24 LB36
## (1,1,0)(1,0,0)[12] 0.07489526 0.2695822
## (1,1,0)(0,0,1)[12] 0.06336111 0.2285484
Os testes de Ljung-box para os lags \(24\) e \(36\) apresentaram p-valores superiores a \(5\%\), o que indica ausênsia de correlação residual, e então que os modelos estão capturando adequadamente a estrutura temporal da série, sem restar um padrão sistemático claro dos erros aparentemente.